package scu.maqiang.homogenization;

import scu.maqiang.fes.*;
import scu.maqiang.fes.FES2T32;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.*;

public class ElasticitySOTS2DT3 {

	public static void main(String[] args) {
		long begin = System.nanoTime();
	    int N = 8;
	    double eps= 1.0 / N;
	    Mesh2T3 cellMesh = new Mesh2T3().square(64, 64);
	    ScalarFunc regionCellFunc = (xy, llable, param) -> {
            return (Math.abs(xy[0] - 0.5) < 0.25 && Math.abs(xy[1] - 0.5) < 0.25)? 1.0:-1.0;
        };
	    cellMesh.setDomainLabel(regionCellFunc, null, 10);
	    FES2T32 fsCell = new FES2T32(cellMesh);
	    double[][] FOCS = new double[3][fsCell.GetNdof()];
		double[][] SOCS = new double[6][fsCell.GetNdof()];
		double[][] D = new double[3][3];
		SOTS2DT362 ss = new SOTS2DT362(fsCell, null);
		
		//ScalarCoefFunc[] gcf = {cf};
		
		double E1 = 1.5e4, E2 = 2.01e6;
		double Nu1 = 0.3, Nu2 = 0.25;
		double[][] D1 = Materials.ElasticConstituteMatrix(E1, Nu1, BVPType.PLANESTRAIN);
		double[][] D2 = Materials.ElasticConstituteMatrix(E2, Nu2, BVPType.PLANESTRAIN);
		MatrixFunc mf = (xyz, label, param) -> label == 10? D1 : D2;
		ss.elasticityCellSolution(mf, null, BVPType.CONSTITUTE_MATRIX_COMMON, FOCS, D, SOCS);
	    
		System.out.println("min: " + MVO.min(FOCS[0]) + "\tmax: " + MVO.max(FOCS[0]));
		System.out.println("min: " + MVO.min(FOCS[1]) + "\tmax: " + MVO.max(FOCS[1]));
		System.out.println("min: " + MVO.min(FOCS[2]) + "\tmax: " + MVO.max(FOCS[2]));
		System.out.println("Homogenized Elasticity: ");
		System.out.println(MVO.toString(D));
		System.out.println("min: " + MVO.min(SOCS[0]) + "\tmax: " + MVO.max(SOCS[0]));
		System.out.println("min: " + MVO.min(SOCS[1]) + "\tmax: " + MVO.max(SOCS[1]));
		System.out.println("min: " + MVO.min(SOCS[2]) + "\tmax: " + MVO.max(SOCS[2]));
		System.out.println("min: " + MVO.min(SOCS[3]) + "\tmax: " + MVO.max(SOCS[3]));
		System.out.println("min: " + MVO.min(SOCS[4]) + "\tmax: " + MVO.max(SOCS[4]));
		System.out.println("min: " + MVO.min(SOCS[5]) + "\tmax: " + MVO.max(SOCS[5]));
		//cellMesh.toTecplot("CellFOCS.dat", FOCS);
		//cellMesh.toTecplot("CellSOCS.dat", SOCS);

		Mesh2T6 homoMesh = new Mesh2T6().square(256, 32, T3Type.Cross2);
		homoMesh.scale(8.0, 1.0);
		FES2T62 fsHomo = new FES2T62(homoMesh);
		D[0][2] = D[1][2] = D[2][0] = D[2][1] = 0.0;
		MatrixFunc homoCoef = (xy, label, param) -> D;
		SRMatrix homoA = new SRMatrix(fsHomo.GetNdof());
		fsHomo.assembleStiff(homoCoef, null, BVPType.CONSTITUTE_MATRIX_COMMON, homoA);
		double[] homoRHS = new double[fsHomo.GetNdof()];
		fsHomo.assembleSource(new double[] {0, -1.0e3}, BVPType.PLANESTRAIN, homoRHS);
		fsHomo.applyBC_MBN(homoA, Direct.All, 4);
		fsHomo.applyBC_MBN(homoRHS, Direct.All, 0.0, 4);
		
//		NewIterSSolver homoSolver = new NewIterSSolver(homoA);
//		double[] homoSol = new double[fsHomo.GetNdof()];
//		homoSolver.PCGSSOR(homoRHS, homoSol, 1.5, 1);

		DirectSSolver homoSolver = new DirectSSolver(homoA);
		double[] homoSol = homoSolver.CholeskySolve(homoRHS, 1);

		double[][] homoSol_uv = new double[2][homoMesh.getNv()];
		fsHomo.extractComponent(homoSol, homoSol_uv);
		
        System.out.println("Homo Solution");
        System.out.println("u0 min: " + MVO.min(homoSol_uv[0]) + "\tmax: " + MVO.max(homoSol_uv[0]));
        System.out.println("v0 min: " + MVO.min(homoSol_uv[1]) + "\tmax: " + MVO.max(homoSol_uv[1]));
		
		Mesh2T3 fineMesh = new Mesh2T3().square(2048, 256);
		fineMesh.scale(8.0, 1.0);
		ScalarFunc regionFineFunc = (xy, rrlabel, param) -> {
            double xCell = xy[0] * N - (int)(xy[0] * N);
            double yCell = xy[1] * N - (int)(xy[1] * N);
            return (Math.abs(xCell - 0.5) < 0.25 && Math.abs(yCell - 0.5) < 0.25)? 1.0:-1.0;
        };
        fineMesh.setDomainLabel(regionFineFunc, null, 10);
        FES2T32 fsFine = new FES2T32(fineMesh);
        SRMatrix fineA = new SRMatrix(fsFine.GetNdof());
        fsFine.assembleElasticStiff(mf, null, BVPType.CONSTITUTE_MATRIX_COMMON, fineA);
        double[] fineRHS = new double[fsFine.GetNdof()];
        fsFine.assembleVolumeForce(new double[] {0, -1.0e3}, BVPType.PLANESTRAIN, fineRHS);
        fsFine.applyBC_MBN(fineA, Direct.All, 4);
        fsFine.applyBC_MBN(fineRHS, Direct.All, 0.0, 4);
        
        //NewIterSSolver fineSolver = new NewIterSSolver(fineA);
        //double[] fineSol = new double[fsFine.GetNdof()];
        //fineSolver.PCGSSOR(fineRHS, fineSol, 1.5, 1);
		DirectSSolver fineSolver = new DirectSSolver(fineA);
		double[] fineSol = fineSolver.CholeskySolve(fineRHS, 1);

        double[][] fineSol_uv = new double[2][fineMesh.getNv()];
        fsFine.extractComponent(fineSol, fineSol_uv);
        
        System.out.println("Fine Solution");
        System.out.println("u min: " + MVO.min(fineSol_uv[0]) + "\tmax: " + MVO.max(fineSol_uv[0]));
        System.out.println("v min: " + MVO.min(fineSol_uv[1]) + "\tmax: " + MVO.max(fineSol_uv[1]));
        
        double[][] HomoFineSol_uv = new double[2][fineMesh.getNv()];
        double[][] FirstSol_uv = new double[2][fineMesh.getNv()];
        double[][] SecondSol_uv = new double[2][fineMesh.getNv()];
        
        ss.assembleElasticitySOTSSolution(fsFine, HomoFineSol_uv, FirstSol_uv, SecondSol_uv,
        		                          fsHomo, homoSol_uv, fsCell, FOCS, SOCS, N, eps);
        System.out.println("Homo Solution");
        System.out.println("u0 min: " + MVO.min(HomoFineSol_uv[0]) + "\tmax: " + MVO.max(HomoFineSol_uv[0]));
        System.out.println("v0 min: " + MVO.min(HomoFineSol_uv[1]) + "\tmax: " + MVO.max(HomoFineSol_uv[1]));
        System.out.println("First Order Approximation: ");
        System.out.println("u1 min: " + MVO.min(FirstSol_uv[0]) + "\tmax: " + MVO.max(FirstSol_uv[0]));
        System.out.println("v1 min: " + MVO.min(FirstSol_uv[1]) + "\tmax: " + MVO.max(FirstSol_uv[1]));
        System.out.println("Second Order Approximation: ");
        System.out.println("u2 min: " + MVO.min(SecondSol_uv[0]) + "\tmax: " + MVO.max(SecondSol_uv[0]));
        System.out.println("v2 min: " + MVO.min(SecondSol_uv[1]) + "\tmax: " + MVO.max(SecondSol_uv[1]));
        
        double[][] wholeSol = new double[8][];
        wholeSol[0] = fineSol_uv[0];
        wholeSol[1] = fineSol_uv[1];
        wholeSol[2] = HomoFineSol_uv[0];
        wholeSol[3] = HomoFineSol_uv[1];
        wholeSol[4] = FirstSol_uv[0];
        wholeSol[5] = FirstSol_uv[1];
        wholeSol[6] = SecondSol_uv[0];
        wholeSol[7] = SecondSol_uv[1];
        fineMesh.toTecplot("ElasticSOTSSolutions.dat", wholeSol);
        
        double[] error0_u = MVO.add(fineSol_uv[0], -1.0, HomoFineSol_uv[0]);
		double[] error1_u = MVO.add(fineSol_uv[0], -1.0, FirstSol_uv[0]);
		double[] error2_u = MVO.add(fineSol_uv[0], -1.0, SecondSol_uv[0]);
		FES2T31 fs1 = new FES2T31(fineMesh);
		double L2Norm = fs1.computeL2Norm(fineSol_uv[0]);
		System.out.println("L2Rrror:");
		System.out.println("u Error 0: " + fs1.computeL2Norm(error0_u) / L2Norm);
		System.out.println("u Error 1: " + fs1.computeL2Norm(error1_u) / L2Norm);
		System.out.println("u Error 2: " + fs1.computeL2Norm(error2_u) / L2Norm);
		double H1Norm = fs1.computeH1SemiNorm(fineSol_uv[0]);
		System.out.println("H1Rrror:");
		System.out.println("u Error 0: " + fs1.computeH1SemiNorm(error0_u) / H1Norm);
		System.out.println("u Error 1: " + fs1.computeH1SemiNorm(error1_u) / H1Norm);
		System.out.println("u Error 2: " + fs1.computeH1SemiNorm(error2_u) / H1Norm);
		
        double[] error0_v = MVO.add(fineSol_uv[1], -1.0, HomoFineSol_uv[1]);
		double[] error1_v = MVO.add(fineSol_uv[1], -1.0, FirstSol_uv[1]);
		double[] error2_v = MVO.add(fineSol_uv[1], -1.0, SecondSol_uv[1]);
		L2Norm = fs1.computeL2Norm(fineSol_uv[1]);
		System.out.println("L2Rrror:");
		System.out.println("v Error 0: " + fs1.computeL2Norm(error0_v) / L2Norm);
		System.out.println("v Error 1: " + fs1.computeL2Norm(error1_v) / L2Norm);
		System.out.println("v Error 2: " + fs1.computeL2Norm(error2_v) / L2Norm);
		H1Norm = fs1.computeH1SemiNorm(fineSol_uv[1]);
		System.out.println("H1Rrror:");
		System.out.println("v Error 0: " + fs1.computeH1SemiNorm(error0_v) / H1Norm);
		System.out.println("v Error 1: " + fs1.computeH1SemiNorm(error1_v) / H1Norm);
		System.out.println("v Error 2: " + fs1.computeH1SemiNorm(error2_v) / H1Norm);
		
		long end = System.nanoTime();
		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");

	}

}
